2  PBMCs: Compile & Format Data

Some longitudinal PBMC timepoints were relabelled after this analysis was conducted. These are the corrected timepoints: P106: W48 -> W49 P109: W48 -> W49 P110: W38 -> W48 P111: W40 -> W41

These timepoints are synonyms: Pre-vaccine <-> Post-nivolumab

2.1 Set up workspace

# Libraries
library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats   1.0.0     ✔ readr     2.1.5
✔ ggplot2   3.5.1     ✔ stringr   1.5.1
✔ lubridate 1.9.4     ✔ tibble    3.2.1
✔ purrr     1.0.4     ✔ tidyr     1.3.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

2.2 Set pseudo count for in vitro samples

Set a pseudo count to allow for in vitro fold change calculations of undetected clones The appropriate value was explored in file “Test_pseudocount_frequency.Rmd” and the slides titled “Testing pseudocounts in calculating fold change”

pseudo_freq <- 0.001

2.3 Open TCR reactivity data

reactivity <- read.csv("/jsimonlab/users/chloetu/melanoma/Melanoma_Eryn/analysis_093024/Eryn_reactivity_results/reformatted_reconstructed_TCR_072324.csv") %>%
  mutate(Beta_clonotype = paste0(TRBV_1, ";", TRBJ_1, ";", CDR3B_1)) %>%
  select(Beta_clonotype, Patient, Reactive) %>%
  filter(Reactive == "Yes")

reactivity %>%
  dplyr::count(Patient)
  Patient  n
1     101  3
2     103  3
3     104 18
4     108  7

2.4 Compile P101’s data

Load longitudinal frequencies

# Control == pretreatment
p101_pretreatment_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/101pre-treament_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

# Post-Nivo == prevax
p101_prevax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/101pre-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

# Post-vaccine
p101_postvax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/101post-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

Aggregate beta chains into one dataframe

p101_betas_pbmc <- p101_pretreatment_beta %>%
  full_join(p101_prevax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p101_pretreatment" = "Freq.x",
                "p101_prevax" = "Freq.y",
                "p101_pretreatment_umi" = "totalUMICount.x",
                "p101_prevax_umi" = "totalUMICount.y") %>%
  full_join(p101_postvax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p101_postvax" = "Freq",
                "p101_postvax_umi" = "totalUMICount") %>%
  select("v_hit", "j_hit", "cdr3", "p101_pretreatment", "p101_pretreatment_umi", 
         "p101_prevax", "p101_prevax_umi", "p101_postvax", "p101_postvax_umi") %>%
  mutate_at(c("p101_pretreatment", "p101_prevax", "p101_postvax"), as.numeric) %>% 
  unite("Beta_clonotype", c(v_hit, j_hit, cdr3), sep = ";")
  

dim(p101_betas_pbmc)
[1] 413697      7

Calculate fold changes across longitudinal timepoints

p101_betas_pbmc <- p101_betas_pbmc %>%
  mutate(postvax_vs_prevax_fc = p101_postvax/p101_prevax,
         prevax_vs_pretreatment_fc = p101_prevax/p101_pretreatment,
         log2fc_postvax_vs_prevax = log2(postvax_vs_prevax_fc),
         log2fc_prevax_vs_pretreatment = log2(prevax_vs_pretreatment_fc))

Load in vitro frequencies

p101DMSO_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/101-DMSO_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p101poolA_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/101-PoolA_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p101poolB_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/101-PoolB_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")

Aggregate beta chains into one dataframe

p101_betas_invitro <- full_join(p101DMSO_beta,p101poolA_beta,by="Beta_clonotype") %>%
    dplyr::rename("p101_DMSO_beta" = Freq.x, 
                  "p101_DMSO_umi" = totalUMICount.x, 
                  "p101_poolA_beta" = Freq.y,
                  "p101_poolA_umi" = totalUMICount.y) %>%
    full_join(p101poolB_beta,by="Beta_clonotype") %>%
    dplyr::rename("p101_poolB_beta" = Freq,
                  "p101_poolB_umi" = totalUMICount)

p101_betas_invitro %>% distinct(Beta_clonotype) %>% dplyr::count()
       n
1 241961

Convert to decimal frequency to %

p101_betas_pbmc <- p101_betas_pbmc %>%
  mutate(p101_pretreatment = p101_pretreatment * 100,
         p101_postvax = p101_postvax * 100,
         p101_prevax = p101_prevax * 100)

p101_betas_invitro <- p101_betas_invitro %>%
  mutate(p101_DMSO_beta = p101_DMSO_beta * 100,
         p101_poolA_beta = p101_poolA_beta * 100,
         p101_poolB_beta = p101_poolB_beta * 100)

Join longitudinal and in vitro frequencies by beta chain

p101_betas <- full_join(p101_betas_pbmc, p101_betas_invitro)
Joining with `by = join_by(Beta_clonotype)`

Calculate fold changes across in vitro samples

p101_betas <- p101_betas %>%
  rowwise() %>%
  mutate(poolA_vs_dmso_fc = sum(c(p101_poolA_beta, pseudo_freq), na.rm = TRUE)/sum(c(p101_DMSO_beta, pseudo_freq), na.rm = TRUE),
         poolB_vs_dmso_fc = sum(c(p101_poolB_beta, pseudo_freq), na.rm = TRUE)/sum(c(p101_DMSO_beta, pseudo_freq), na.rm = TRUE)) %>%
         ungroup() %>%
  mutate(log2fc_poolA_vs_dmso = log2(poolA_vs_dmso_fc),
         log2fc_poolB_vs_dmso = log2(poolB_vs_dmso_fc),
         max_log2fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ log2fc_poolA_vs_dmso,
                                             log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ log2fc_poolB_vs_dmso),
         max_fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ poolA_vs_dmso_fc,
                                         log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ poolB_vs_dmso_fc))

Attach TCR reactivities to dataframe via the beta clonotype

p101_reactive <- reactivity %>%
  filter(Patient == "101")

p101_betas <- p101_betas %>%
  mutate(Reactive = case_when(Beta_clonotype %in% p101_reactive$Beta_clonotype ~ TRUE,
                              !(Beta_clonotype %in% p101_reactive$Beta_clonotype) ~ FALSE))

Save data

write.csv(p101_betas, "p101_betas_merged_Part1.csv", row.names = FALSE)

2.5 Compile P103’s data

Load longitudinal frequencies

p103_pretreatment_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/103pre-treatment_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p103_prevax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/103pre-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p103_postvax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/103post-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p103_w48_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/103wk48_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") 

p103_w72_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/103wk72_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

Aggregate beta chains into one dataframe

p103_betas_pbmc <- p103_pretreatment_beta %>%
  full_join(p103_prevax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p103_pretreatment" = "Freq.x",
                "p103_prevax" = "Freq.y",
                "p103_pretreatment_umi" = "totalUMICount.x",
                "p103_prevax_umi" = "totalUMICount.y") %>%
  full_join(p103_postvax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p103_postvax" = "Freq",
                "p103_postvax_umi" = "totalUMICount") %>% 
  full_join(p103_w48_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p103_w48" = "Freq",
                "p103_w48_umi" = "totalUMICount") %>% 
  full_join(p103_w72_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p103_w72" = "Freq",
                "p103_w72_umi" = "totalUMICount") %>% 
  select("v_hit", "j_hit", "cdr3", "p103_pretreatment", "p103_pretreatment_umi", "p103_prevax", "p103_prevax_umi", "p103_postvax", "p103_postvax_umi", "p103_w48", "p103_w48_umi", "p103_w72", "p103_w72_umi") %>%
  mutate_at(c("p103_pretreatment", "p103_pretreatment_umi", "p103_prevax", "p103_prevax_umi", "p103_postvax", "p103_postvax_umi", "p103_w48", "p103_w48_umi", "p103_w72", "p103_w72_umi"), as.numeric) %>% 
  unite("Beta_clonotype", c(v_hit, j_hit, cdr3), sep = ";")
  

dim(p103_betas_pbmc)
[1] 324551     11

Calculate fold changes across longitudinal timepoints

p103_betas_pbmc <- p103_betas_pbmc %>%
  mutate(postvax_vs_prevax_fc = p103_postvax/p103_prevax,
         prevax_vs_pretreatment_fc = p103_prevax/p103_pretreatment,
         log2fc_postvax_vs_prevax = log2(postvax_vs_prevax_fc),
         log2fc_prevax_vs_pretreatment = log2(prevax_vs_pretreatment_fc))

Load in vitro frequencies

p103DMSO_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/103-DMSO_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p103poolA_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/103-PoolA_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p103poolB_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/103-PoolB_TRB_report.tsv", sep = "\t", row.names = 1) %>%
select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")

Aggregate beta chains into one dataframe

p103_betas_invitro <- full_join(p103DMSO_beta,p103poolA_beta,by="Beta_clonotype") %>%
    dplyr::rename("p103_DMSO_beta" = Freq.x, 
                  "p103_DMSO_umi" = totalUMICount.x, 
                  "p103_poolA_beta" = Freq.y,
                  "p103_poolA_umi" = totalUMICount.y) %>%
    full_join(p103poolB_beta,by="Beta_clonotype") %>%
    dplyr::rename("p103_poolB_beta" = Freq,
                  "p103_poolB_umi" = totalUMICount)

p103_betas_invitro %>% distinct(Beta_clonotype) %>% dplyr::count()
       n
1 128158

Convert to decimal frequency to %

p103_betas_pbmc <- p103_betas_pbmc %>%
  mutate(p103_pretreatment = p103_pretreatment * 100,
         p103_postvax = p103_postvax * 100,
         p103_prevax = p103_prevax * 100,
         p103_w48 = p103_w48 * 100,
         p103_w72 = p103_w72 * 100)

p103_betas_invitro <- p103_betas_invitro %>%
  mutate(p103_DMSO_beta = p103_DMSO_beta * 100,
         p103_poolA_beta = p103_poolA_beta * 100,
         p103_poolB_beta = p103_poolB_beta * 100)

Join longitudinal and in vitro frequencies by beta chain

p103_betas <- full_join(p103_betas_pbmc, p103_betas_invitro)
Joining with `by = join_by(Beta_clonotype)`

Calculate fold changes across in vitro samples

p103_betas <- p103_betas %>%
  rowwise() %>%
  mutate(poolA_vs_dmso_fc = sum(c(p103_poolA_beta, pseudo_freq), na.rm = TRUE)/sum(c(p103_DMSO_beta, pseudo_freq), na.rm = TRUE),
         poolB_vs_dmso_fc = sum(c(p103_poolB_beta, pseudo_freq), na.rm = TRUE)/sum(c(p103_DMSO_beta, pseudo_freq), na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(log2fc_poolA_vs_dmso = log2(poolA_vs_dmso_fc),
         log2fc_poolB_vs_dmso = log2(poolB_vs_dmso_fc),
         max_log2fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ log2fc_poolA_vs_dmso,
                                             log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ log2fc_poolB_vs_dmso),
         max_fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ poolA_vs_dmso_fc,
                                         log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ poolB_vs_dmso_fc))

Attach TCR reactivities to dataframe via the beta clonotype

p103_reactive <- reactivity %>%
  filter(Patient == "103")

p103_betas <- p103_betas %>%
  mutate(Reactive = case_when(Beta_clonotype %in% p103_reactive$Beta_clonotype ~ TRUE,
                              !(Beta_clonotype %in% p103_reactive$Beta_clonotype) ~ FALSE))

Save data

write.csv(p103_betas, "p103_betas_merged_Part1.csv", row.names = FALSE)

2.6 Compile P104’s data

Load longitudinal frequencies

p104_pretreatment_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/104pre-treament_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p104_prevax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/104pre-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p104_postvax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/104post-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p104_w48_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/104wk48_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

Aggregate beta chains into one dataframe

p104_betas_pbmc <- p104_pretreatment_beta %>%
  full_join(p104_prevax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p104_pretreatment" = "Freq.x",
                "p104_prevax" = "Freq.y",
                "p104_pretreatment_umi" = "totalUMICount.x",
                "p104_prevax_umi" = "totalUMICount.y") %>%
  full_join(p104_postvax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p104_postvax" = "Freq",
                "p104_postvax_umi" = "totalUMICount") %>% 
  full_join(p104_w48_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p104_w48" = "Freq",
                "p104_w48_umi" = "totalUMICount") %>% 
  select("v_hit", "j_hit", "cdr3", "p104_pretreatment", "p104_pretreatment_umi", "p104_prevax", "p104_prevax_umi", "p104_postvax", "p104_postvax_umi", "p104_w48", "p104_w48_umi") %>%
  mutate_at(c("p104_pretreatment", "p104_pretreatment_umi", "p104_prevax", "p104_prevax_umi", "p104_postvax", "p104_postvax_umi", "p104_w48", "p104_w48_umi"), as.numeric) %>% 
  unite("Beta_clonotype", c(v_hit, j_hit, cdr3), sep = ";")

dim(p104_betas_pbmc)
[1] 747790      9

Calculate fold changes across longitudinal timepoints

p104_betas_pbmc <- p104_betas_pbmc %>%
  mutate(postvax_vs_prevax_fc = p104_postvax/p104_prevax,
         prevax_vs_pretreatment_fc = p104_prevax/p104_pretreatment,
         log2fc_postvax_vs_prevax = log2(postvax_vs_prevax_fc),
         log2fc_prevax_vs_pretreatment = log2(prevax_vs_pretreatment_fc))

Load in vitro frequencies

p104DMSO_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/104-DMSO_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p104poolA_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/104-PoolA_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p104poolB_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/104-PoolB_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")

Aggregate beta chains into one dataframe

p104_betas_invitro <- full_join(p104DMSO_beta,p104poolA_beta,by="Beta_clonotype") %>%
    dplyr::rename("p104_DMSO_beta" = Freq.x, 
                  "p104_DMSO_umi" = totalUMICount.x, 
                  "p104_poolA_beta" = Freq.y,
                  "p104_poolA_umi" = totalUMICount.y) %>%
    full_join(p104poolB_beta,by="Beta_clonotype") %>%
    dplyr::rename("p104_poolB_beta" = Freq,
                  "p104_poolB_umi" = totalUMICount)

p104_betas_invitro %>% distinct(Beta_clonotype) %>% dplyr::count()
       n
1 440184

Convert to decimal frequency to %

p104_betas_pbmc <- p104_betas_pbmc %>%
  mutate(p104_pretreatment = p104_pretreatment * 100,
         p104_postvax = p104_postvax * 100,
         p104_prevax = p104_prevax * 100,
         p104_w48 = p104_w48 * 100)

p104_betas_invitro <- p104_betas_invitro %>%
  mutate(p104_DMSO_beta = p104_DMSO_beta * 100,
         p104_poolA_beta = p104_poolA_beta * 100,
         p104_poolB_beta = p104_poolB_beta * 100)

Join longitudinal and in vitro frequencies by beta chain

p104_betas <- full_join(p104_betas_pbmc, p104_betas_invitro)
Joining with `by = join_by(Beta_clonotype)`

Calculate fold changes across in vitro samples

p104_betas <- p104_betas %>%
  rowwise() %>%
  mutate(poolA_vs_dmso_fc = sum(c(p104_poolA_beta, pseudo_freq), na.rm = TRUE)/sum(c(p104_DMSO_beta, pseudo_freq), na.rm = TRUE),
         poolB_vs_dmso_fc = sum(c(p104_poolB_beta, pseudo_freq), na.rm = TRUE)/sum(c(p104_DMSO_beta, pseudo_freq), na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(log2fc_poolA_vs_dmso = log2(poolA_vs_dmso_fc),
         log2fc_poolB_vs_dmso = log2(poolB_vs_dmso_fc),
         max_log2fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ log2fc_poolA_vs_dmso,
                                             log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ log2fc_poolB_vs_dmso),
         max_fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ poolA_vs_dmso_fc,
                                         log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ poolB_vs_dmso_fc))

Attach TCR reactivities to dataframe via the beta clonotype

p104_reactive <- reactivity %>%
  filter(Patient == "104")

p104_betas <- p104_betas %>%
  mutate(Reactive = case_when(Beta_clonotype %in% p104_reactive$Beta_clonotype ~ TRUE,
                              !(Beta_clonotype %in% p104_reactive$Beta_clonotype) ~ FALSE))

Save data

write.csv(p104_betas, "p104_betas_merged_Part1.csv", row.names = FALSE)

2.7 Compile P105’s data

Load longitudinal frequencies

p105_pretreatment_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/105pre-treament_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p105_prevax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/105pre-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p105_postvax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/105post-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p105_W48_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/105wk48_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

Aggregate beta chains into one dataframe

p105_betas_pbmc <- p105_pretreatment_beta %>%
  full_join(p105_prevax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p105_pretreatment" = "Freq.x",
                "p105_prevax" = "Freq.y",
                "p105_pretreatment_umi" = "totalUMICount.x",
                "p105_prevax_umi" = "totalUMICount.y") %>%
  full_join(p105_postvax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p105_postvax" = "Freq",
                "p105_postvax_umi" = "totalUMICount") %>% 
  full_join(p105_W48_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p105_w48" = "Freq",
                "p105_w48_umi" = "totalUMICount") %>%
  select("v_hit", "j_hit", "cdr3", "p105_pretreatment", "p105_pretreatment_umi", "p105_prevax", "p105_prevax_umi", "p105_postvax", "p105_postvax_umi", "p105_w48", "p105_w48_umi") %>%
  mutate_at(c("p105_pretreatment", "p105_pretreatment_umi", "p105_prevax", "p105_prevax_umi", "p105_postvax", "p105_postvax_umi", "p105_w48", "p105_w48_umi"), as.numeric) %>% 
  unite("Beta_clonotype", c(v_hit, j_hit, cdr3), sep = ";")

dim(p105_betas_pbmc)
[1] 465983      9

Calculate fold changes across longitudinal timepoints

p105_betas_pbmc <- p105_betas_pbmc %>%
  mutate(postvax_vs_prevax_fc = p105_postvax/p105_prevax,
         prevax_vs_pretreatment_fc = p105_prevax/p105_pretreatment,
         log2fc_postvax_vs_prevax = log2(postvax_vs_prevax_fc),
         log2fc_prevax_vs_pretreatment = log2(prevax_vs_pretreatment_fc))

Load in vitro frequencies

p105DMSO_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/105-DMSO_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p105poolA_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/105-PoolA_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p105poolB_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/105-PoolB_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")

Aggregate beta chains into one dataframe

p105_betas_invitro <- full_join(p105DMSO_beta,p105poolA_beta,by="Beta_clonotype") %>%
    dplyr::rename("p105_DMSO_beta" = Freq.x, 
                  "p105_DMSO_umi" = totalUMICount.x, 
                  "p105_poolA_beta" = Freq.y,
                  "p105_poolA_umi" = totalUMICount.y) %>%
    full_join(p105poolB_beta,by="Beta_clonotype") %>%
    dplyr::rename("p105_poolB_beta" = Freq,
                  "p105_poolB_umi" = totalUMICount)

p105_betas_invitro %>% distinct(Beta_clonotype) %>% dplyr::count()
       n
1 182772

Convert to decimal frequency to %

p105_betas_pbmc <- p105_betas_pbmc %>%
  mutate(p105_pretreatment = p105_pretreatment * 100,
         p105_postvax = p105_postvax * 100,
         p105_prevax = p105_prevax * 100,
         p105_w48 = p105_w48 * 100)

p105_betas_invitro <- p105_betas_invitro %>%
  mutate(p105_DMSO_beta = p105_DMSO_beta * 100,
         p105_poolA_beta = p105_poolA_beta * 100,
         p105_poolB_beta = p105_poolB_beta * 100)

Join longitudinal and in vitro frequencies by beta chain

p105_betas <- full_join(p105_betas_pbmc, p105_betas_invitro)
Joining with `by = join_by(Beta_clonotype)`

Calculate fold changes across in vitro samples

p105_betas <- p105_betas %>%
  rowwise() %>%
  mutate(poolA_vs_dmso_fc = sum(c(p105_poolA_beta, pseudo_freq), na.rm = TRUE)/sum(c(p105_DMSO_beta, pseudo_freq), na.rm = TRUE),
         poolB_vs_dmso_fc = sum(c(p105_poolB_beta, pseudo_freq), na.rm = TRUE)/sum(c(p105_DMSO_beta, pseudo_freq), na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(log2fc_poolA_vs_dmso = log2(poolA_vs_dmso_fc),
         log2fc_poolB_vs_dmso = log2(poolB_vs_dmso_fc),
         max_log2fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ log2fc_poolA_vs_dmso,
                                             log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ log2fc_poolB_vs_dmso),
         max_fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ poolA_vs_dmso_fc,
                                         log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ poolB_vs_dmso_fc))

Save data

write.csv(p105_betas, "p105_betas_merged_Part1.csv", row.names = FALSE)

2.8 Compile P106’s data

Load longitudinal frequencies

p106_pretreatment_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/106pre-treament_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p106_prevax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/106pre-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p106_postvax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/106post-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p106_w48_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/106wk48_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

Aggregate beta chains into one dataframe

p106_betas_pbmc <- p106_pretreatment_beta %>%
  full_join(p106_prevax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p106_pretreatment" = "Freq.x",
                "p106_prevax" = "Freq.y",
                "p106_pretreatment_umi" = "totalUMICount.x",
                "p106_prevax_umi" = "totalUMICount.y") %>%
  full_join(p106_postvax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p106_postvax" = "Freq",
                "p106_postvax_umi" = "totalUMICount") %>% 
  full_join(p106_w48_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p106_w48" = "Freq",
                "p106_w48_umi" = "totalUMICount") %>% 
  select("v_hit", "j_hit", "cdr3", "p106_pretreatment", "p106_pretreatment_umi", "p106_prevax", "p106_prevax_umi", "p106_postvax", "p106_postvax_umi", "p106_w48", "p106_w48_umi") %>%
  mutate_at(c("p106_pretreatment", "p106_pretreatment_umi", "p106_prevax", "p106_prevax_umi", "p106_postvax", "p106_postvax_umi", "p106_w48", "p106_w48_umi"), as.numeric) %>% 
  unite("Beta_clonotype", c(v_hit, j_hit, cdr3), sep = ";")
  

dim(p106_betas_pbmc)
[1] 662674      9

Calculate fold changes across longitudinal timepoints

p106_betas_pbmc <- p106_betas_pbmc %>%
  mutate(postvax_vs_prevax_fc = p106_postvax/p106_prevax,
         prevax_vs_pretreatment_fc = p106_prevax/p106_pretreatment,
         log2fc_postvax_vs_prevax = log2(postvax_vs_prevax_fc),
         log2fc_prevax_vs_pretreatment = log2(prevax_vs_pretreatment_fc))

Load in vitro frequencies

p106DMSO_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/106-DMSO_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p106poolA_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/106-PoolA_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p106poolB_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/106-PoolB_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")

Aggregate beta chains into one dataframe

p106_betas_invitro <- full_join(p106DMSO_beta,p106poolA_beta,by="Beta_clonotype") %>%
    dplyr::rename("p106_DMSO_beta" = Freq.x, 
                  "p106_DMSO_umi" = totalUMICount.x, 
                  "p106_poolA_beta" = Freq.y,
                  "p106_poolA_umi" = totalUMICount.y) %>%
    full_join(p106poolB_beta,by="Beta_clonotype") %>%
    dplyr::rename("p106_poolB_beta" = Freq,
                  "p106_poolB_umi" = totalUMICount)

p106_betas_invitro %>% distinct(Beta_clonotype) %>% dplyr::count()
       n
1 341472

Convert to decimal frequency to %

p106_betas_pbmc <- p106_betas_pbmc %>%
  mutate(p106_pretreatment = p106_pretreatment * 100,
         p106_postvax = p106_postvax * 100,
         p106_prevax = p106_prevax * 100,
         p106_w48 = p106_w48 * 100)

p106_betas_invitro <- p106_betas_invitro %>%
  mutate(p106_DMSO_beta = p106_DMSO_beta * 100,
         p106_poolA_beta = p106_poolA_beta * 100,
         p106_poolB_beta = p106_poolB_beta * 100)

Join longitudinal and in vitro frequencies by beta chain

p106_betas <- full_join(p106_betas_pbmc, p106_betas_invitro)
Joining with `by = join_by(Beta_clonotype)`

Calculate fold changes across in vitro samples

p106_betas <- p106_betas %>%
  rowwise() %>%
  mutate(poolA_vs_dmso_fc = sum(c(p106_poolA_beta, pseudo_freq), na.rm = TRUE)/sum(c(p106_DMSO_beta, pseudo_freq), na.rm = TRUE),
         poolB_vs_dmso_fc = sum(c(p106_poolB_beta, pseudo_freq), na.rm = TRUE)/sum(c(p106_DMSO_beta, pseudo_freq), na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(log2fc_poolA_vs_dmso = log2(poolA_vs_dmso_fc),
         log2fc_poolB_vs_dmso = log2(poolB_vs_dmso_fc),
         max_log2fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ log2fc_poolA_vs_dmso,
                                             log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ log2fc_poolB_vs_dmso),
         max_fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ poolA_vs_dmso_fc,
                                         log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ poolB_vs_dmso_fc))

Save data

write.csv(p106_betas, "p106_betas_merged_Part1.csv", row.names = FALSE)

2.9 Compile P108’s data

Load longitudinal frequencies

p108_pretreatment_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/108pre-treatment_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p108_prevax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/108pre-vaccine_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p108_postvax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/108post-vaccine_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p108_w32_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/108-Wk32d1_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

Aggregate beta chains into one dataframe

p108_betas_pbmc <- p108_pretreatment_beta %>%
  full_join(p108_prevax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p108_pretreatment" = "Freq.x",
                "p108_prevax" = "Freq.y",
                "p108_pretreatment_umi" = "totalUMICount.x",
                "p108_prevax_umi" = "totalUMICount.y"
                ) %>%
  full_join(p108_postvax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p108_postvax" = "Freq",
                "p108_postvax_umi" = "totalUMICount") %>% 
  full_join(p108_w32_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p108_w32" = "Freq",
                "p108_w32_umi" = "totalUMICount") %>% 
  select("v_hit", "j_hit", "cdr3", "p108_pretreatment", "p108_pretreatment_umi", "p108_prevax", "p108_prevax_umi", "p108_postvax", "p108_postvax_umi", "p108_w32", "p108_w32_umi") %>%
  mutate_at(c("p108_pretreatment", "p108_pretreatment_umi", "p108_prevax", "p108_prevax_umi", "p108_postvax", "p108_postvax_umi", "p108_w32", "p108_w32_umi"), as.numeric) %>% 
  unite("Beta_clonotype", c(v_hit, j_hit, cdr3), sep = ";")

  
dim(p108_betas_pbmc)
[1] 674885      9

Calculate fold changes across longitudinal timepoints

p108_betas_pbmc <- p108_betas_pbmc %>%
  mutate(postvax_vs_prevax_fc = p108_postvax/p108_prevax,
         prevax_vs_pretreatment_fc = p108_prevax/p108_pretreatment,
         log2fc_postvax_vs_prevax = log2(postvax_vs_prevax_fc),
         log2fc_prevax_vs_pretreatment = log2(prevax_vs_pretreatment_fc))

Load in vitro frequencies

p108DMSO_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/108-DMSO_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p108poolA_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/108-PoolA-B_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p108poolB_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/108-PoolB_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")

Aggregate beta chains into one dataframe

p108_betas_invitro <- full_join(p108DMSO_beta,p108poolA_beta,by="Beta_clonotype") %>%
    dplyr::rename("p108_DMSO_beta" = Freq.x,
                  "p108_DMSO_umi" = totalUMICount.x, 
                  "p108_poolA_beta" = Freq.y,
                  "p108_poolA_umi" = totalUMICount.y) %>%
    full_join(p108poolB_beta,by="Beta_clonotype") %>%
    dplyr::rename("p108_poolB_beta" = Freq,
                  "p106_poolB_umi" = totalUMICount)

p108_betas_invitro %>% distinct(Beta_clonotype) %>% dplyr::count()
       n
1 315151

Convert to decimal frequency to %

p108_betas_pbmc <- p108_betas_pbmc %>%
  mutate(p108_pretreatment = p108_pretreatment * 100,
         p108_postvax = p108_postvax * 100,
         p108_prevax = p108_prevax * 100,
         p108_w32 = p108_w32 * 100)

p108_betas_invitro <- p108_betas_invitro %>%
  mutate(p108_DMSO_beta = p108_DMSO_beta * 100,
         p108_poolA_beta = p108_poolA_beta * 100,
         p108_poolB_beta = p108_poolB_beta * 100)

Join longitudinal and in vitro frequencies by beta chain

p108_betas <- full_join(p108_betas_pbmc, p108_betas_invitro)
Joining with `by = join_by(Beta_clonotype)`

Calculate fold changes across in vitro samples

p108_betas <- p108_betas %>%
  rowwise() %>%
  mutate(poolA_vs_dmso_fc = sum(c(p108_poolA_beta, pseudo_freq), na.rm = TRUE)/sum(c(p108_DMSO_beta, pseudo_freq), na.rm = TRUE),
         poolB_vs_dmso_fc = sum(c(p108_poolB_beta, pseudo_freq), na.rm = TRUE)/sum(c(p108_DMSO_beta, pseudo_freq), na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(log2fc_poolA_vs_dmso = log2(poolA_vs_dmso_fc),
         log2fc_poolB_vs_dmso = log2(poolB_vs_dmso_fc),
         max_log2fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ log2fc_poolA_vs_dmso,
                                             log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ log2fc_poolB_vs_dmso),
         max_fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ poolA_vs_dmso_fc,
                                         log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ poolB_vs_dmso_fc))

Attach TCR reactivities to dataframe via the beta clonotype

p108_reactive <- reactivity %>%
  filter(Patient == "108")

p108_betas <- p108_betas %>%
  mutate(Reactive = case_when(Beta_clonotype %in% p108_reactive$Beta_clonotype ~ TRUE,
                              !(Beta_clonotype %in% p108_reactive$Beta_clonotype) ~ FALSE))

Save data

write.csv(p108_betas, "p108_betas_merged_Part1.csv", row.names = FALSE)

2.10 Compile P109’s data

Load longitudinal frequencies

p109_pretreatment_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/109pre-treament_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p109_prevax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/109pre-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p109_postvax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/109post-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p109_w48_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/109wk48_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

Aggregate beta chains into one dataframe

p109_betas_pbmc <- p109_pretreatment_beta %>%
  full_join(p109_prevax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p109_pretreatment" = "Freq.x",
                "p109_prevax" = "Freq.y",
                "p109_pretreatment_umi" = "totalUMICount.x",
                "p109_prevax_umi" = "totalUMICount.y") %>%
  full_join(p109_postvax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p109_postvax" = "Freq",
                "p109_postvax_umi" = "totalUMICount") %>% 
  full_join(p109_w48_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p109_w48" = "Freq",
                "p109_w48_umi" = "totalUMICount") %>% 
  select("v_hit", "j_hit", "cdr3", "p109_pretreatment", "p109_pretreatment_umi", "p109_prevax", "p109_prevax_umi", "p109_postvax", "p109_postvax_umi", "p109_w48", "p109_w48_umi") %>%
  mutate_at(c("p109_pretreatment", "p109_pretreatment_umi", "p109_prevax", "p109_prevax_umi", "p109_postvax", "p109_postvax_umi", "p109_w48", "p109_w48_umi"), as.numeric) %>% 
  unite("Beta_clonotype", c(v_hit, j_hit, cdr3), sep = ";")

dim(p109_betas_pbmc)
[1] 680233      9

Calculate fold changes across longitudinal timepoints

p109_betas_pbmc <- p109_betas_pbmc %>%
  mutate(postvax_vs_prevax_fc = p109_postvax/p109_prevax,
         prevax_vs_pretreatment_fc = p109_prevax/p109_pretreatment,
         log2fc_postvax_vs_prevax = log2(postvax_vs_prevax_fc),
         log2fc_prevax_vs_pretreatment = log2(prevax_vs_pretreatment_fc))

Load in vitro frequencies

p109DMSO_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/18279-109DMSO_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p109poolA_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/18279-109PoolA_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p109poolB_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/18279-109PoolB_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")

Aggregate beta chains into one dataframe

p109_betas_invitro <- full_join(p109DMSO_beta,p109poolA_beta,by="Beta_clonotype") %>%
    dplyr::rename("p109_DMSO_beta" = Freq.x, 
                  "p109_poolA_beta" = Freq.y,
                  "p109_DMSO_umi" = totalUMICount.x,
                  "p109_poolA_umi" = totalUMICount.y) %>%
    full_join(p109poolB_beta,by="Beta_clonotype") %>%
    dplyr::rename("p109_poolB_beta" = Freq,
                  "p109_poolB_umi" = totalUMICount)

p109_betas_invitro %>% distinct(Beta_clonotype) %>% dplyr::count()
        n
1 1078940

Convert to decimal frequency to %

p109_betas_pbmc <- p109_betas_pbmc %>%
  mutate(p109_pretreatment = p109_pretreatment * 100,
         p109_postvax = p109_postvax * 100,
         p109_prevax = p109_prevax * 100,
         p109_w48 = p109_w48 * 100)

p109_betas_invitro <- p109_betas_invitro %>%
  mutate(p109_DMSO_beta = p109_DMSO_beta * 100,
         p109_poolA_beta = p109_poolA_beta * 100,
         p109_poolB_beta = p109_poolB_beta * 100)

Join longitudinal and in vitro frequencies by beta chain

p109_betas <- full_join(p109_betas_pbmc, p109_betas_invitro)
Joining with `by = join_by(Beta_clonotype)`

Calculate fold changes across in vitro samples

p109_betas <- p109_betas %>%
  rowwise() %>%
  mutate(poolA_vs_dmso_fc = sum(c(p109_poolA_beta, pseudo_freq), na.rm = TRUE)/sum(c(p109_DMSO_beta, pseudo_freq), na.rm = TRUE),
         poolB_vs_dmso_fc = sum(c(p109_poolB_beta, pseudo_freq), na.rm = TRUE)/sum(c(p109_DMSO_beta, pseudo_freq), na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(log2fc_poolA_vs_dmso = log2(poolA_vs_dmso_fc),
         log2fc_poolB_vs_dmso = log2(poolB_vs_dmso_fc),
         max_log2fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ log2fc_poolA_vs_dmso,
                                             log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ log2fc_poolB_vs_dmso),
         max_fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ poolA_vs_dmso_fc,
                                         log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ poolB_vs_dmso_fc))

Save data

write.csv(p109_betas, "p109_betas_merged_Part1.csv", row.names = FALSE)

2.11 Compile P110’s data

Load longitudinal frequencies

p110_pretreatment_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/110pre-treat_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p110_prevax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/110pre-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p110_postvax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/110post-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p110_w38_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/110wk38_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

Aggregate beta chains into one dataframe

p110_betas_pbmc <- p110_pretreatment_beta %>%
  full_join(p110_prevax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p110_pretreatment" = "Freq.x",
                "p110_prevax" = "Freq.y",
                "p110_pretreatment_umi" = "totalUMICount.x",
                "p110_prevax_umi" = "totalUMICount.y") %>%
  full_join(p110_postvax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p110_postvax" = "Freq",
                "p110_postvax_umi" = "totalUMICount") %>% 
  full_join(p110_w38_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p110_w38" = "Freq",
                "p110_w38_umi" = "totalUMICount") %>% 
  select("v_hit", "j_hit", "cdr3", "p110_pretreatment", "p110_pretreatment_umi", "p110_prevax", "p110_prevax_umi", "p110_postvax", "p110_postvax_umi", "p110_w38", "p110_w38_umi") %>%
  mutate_at(c("p110_pretreatment", "p110_pretreatment_umi", "p110_prevax", "p110_prevax_umi", "p110_postvax", "p110_postvax_umi", "p110_w38", "p110_w38_umi"), as.numeric) %>% 
  unite("Beta_clonotype", c(v_hit, j_hit, cdr3), sep = ";")

dim(p110_betas_pbmc)
[1] 744561      9

Calculate fold changes across longitudinal timepoints

p110_betas_pbmc <- p110_betas_pbmc %>%
  mutate(postvax_vs_prevax_fc = p110_postvax/p110_prevax,
         prevax_vs_pretreatment_fc = p110_prevax/p110_pretreatment,
         log2fc_postvax_vs_prevax = log2(postvax_vs_prevax_fc),
         log2fc_prevax_vs_pretreatment = log2(prevax_vs_pretreatment_fc))

Load in vitro frequencies

p110DMSO_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/110DMSO_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p110poolA_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/110PoolA_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p110poolB_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/110PoolB_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")

Aggregate beta chains into one dataframe

p110_betas_invitro <- full_join(p110DMSO_beta,p110poolA_beta,by="Beta_clonotype") %>%
    dplyr::rename("p110_DMSO_beta" = Freq.x, 
                  "p110_poolA_beta" = Freq.y,
                  "p110_DMSO_umi" = totalUMICount.x,
                  "p110_poolA_umi" = totalUMICount.y) %>%
    full_join(p110poolB_beta,by="Beta_clonotype") %>%
    dplyr::rename("p110_poolB_beta" = Freq,
                  "p110_poolB_umi" = totalUMICount)

p110_betas_invitro %>% distinct(Beta_clonotype) %>% dplyr::count()
       n
1 416916

Convert to decimal frequency to %

p110_betas_pbmc <- p110_betas_pbmc %>%
  mutate(p110_pretreatment = p110_pretreatment * 100,
         p110_postvax = p110_postvax * 100,
         p110_prevax = p110_prevax * 100,
         p110_w38 = p110_w38 * 100)

p110_betas_invitro <- p110_betas_invitro %>%
  mutate(p110_DMSO_beta = p110_DMSO_beta * 100,
         p110_poolA_beta = p110_poolA_beta * 100,
         p110_poolB_beta = p110_poolB_beta * 100)

Join longitudinal and in vitro frequencies by beta chain

p110_betas <- full_join(p110_betas_pbmc, p110_betas_invitro)
Joining with `by = join_by(Beta_clonotype)`

Calculate fold changes across in vitro samples

p110_betas <- p110_betas %>%
  rowwise() %>%
  mutate(poolA_vs_dmso_fc = sum(c(p110_poolA_beta, pseudo_freq), na.rm = TRUE)/sum(c(p110_DMSO_beta, pseudo_freq), na.rm = TRUE),
         poolB_vs_dmso_fc = sum(c(p110_poolB_beta, pseudo_freq), na.rm = TRUE)/sum(c(p110_DMSO_beta, pseudo_freq), na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(log2fc_poolA_vs_dmso = log2(poolA_vs_dmso_fc),
         log2fc_poolB_vs_dmso = log2(poolB_vs_dmso_fc),
         max_log2fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ log2fc_poolA_vs_dmso,
                                             log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ log2fc_poolB_vs_dmso),
         max_fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ poolA_vs_dmso_fc,
                                         log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ poolB_vs_dmso_fc))

Save data

write.csv(p110_betas, "p110_betas_merged_Part1.csv", row.names = FALSE)

2.12 Compile P111’s data

Load longitudinal frequencies

p111_pretreatment_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/111pre-treat_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p111_prevax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/111pre-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p111_postvax_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/111post-vax_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

p111_w40_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/111wk40_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene")

Aggregate beta chains into one dataframe

p111_betas_pbmc <- p111_pretreatment_beta %>%
  full_join(p111_prevax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p111_pretreatment" = "Freq.x",
                "p111_prevax" = "Freq.y",
                "p111_pretreatment_umi" = "totalUMICount.x",
                "p111_prevax_umi" = "totalUMICount.y") %>%
  full_join(p111_postvax_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p111_postvax" = "Freq",
                "p111_postvax_umi" = "totalUMICount") %>% 
  full_join(p111_w40_beta, by = c("v_hit", "j_hit", "cdr3")) %>%
  dplyr::rename("p111_w40" = "Freq",
                "p111_w40_umi" = "totalUMICount") %>% 
  select("v_hit", "j_hit", "cdr3", "p111_pretreatment", "p111_pretreatment_umi", "p111_prevax", "p111_prevax_umi", "p111_postvax", "p111_postvax_umi", "p111_w40", "p111_w40_umi") %>%
  mutate_at(c("p111_pretreatment", "p111_pretreatment_umi", "p111_prevax", "p111_prevax_umi", "p111_postvax", "p111_postvax_umi", "p111_w40", "p111_w40_umi"), as.numeric) %>% 
  unite("Beta_clonotype", c(v_hit, j_hit, cdr3), sep = ";")
  
dim(p111_betas_pbmc)
[1] 607386      9

Calculate fold changes across longitudinal timepoints

p111_betas_pbmc <- p111_betas_pbmc %>%
  mutate(postvax_vs_prevax_fc = p111_postvax/p111_prevax,
         prevax_vs_pretreatment_fc = p111_prevax/p111_pretreatment,
         log2fc_postvax_vs_prevax = log2(postvax_vs_prevax_fc),
         log2fc_prevax_vs_pretreatment = log2(prevax_vs_pretreatment_fc))

Load in vitro frequencies

p111DMSO_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/111-DMSO_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p111poolA_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/111-PoolA_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")
p111poolB_beta <- read.csv("/jsimonlab/projects/Wu/Melanoma_bulkTCR_Eryn/mixcr_processed_031324/Eryn_bulkrhTCR/111-PoolB_TRB_report.tsv", sep = "\t", row.names = 1) %>%
  select(aaSeqCDR3, vGene, jGene, totalUMICount, Freq) %>%
  dplyr::rename("cdr3" = "aaSeqCDR3", "v_hit" = "vGene", "j_hit" = "jGene") %>%
  unite("Beta_clonotype", c(v_hit,j_hit,cdr3),sep=";")

Aggregate beta chains into one dataframe

p111_betas_invitro <- full_join(p111DMSO_beta,p111poolA_beta,by="Beta_clonotype") %>%
    dplyr::rename("p111_DMSO_beta" = Freq.x,
                  "p111_poolA_beta" = Freq.y,
                  "p111_DMSO_umi" = totalUMICount.x,
                  "p111_poolA_umi" = totalUMICount.y) %>%
    full_join(p111poolB_beta,by="Beta_clonotype") %>%
    dplyr::rename("p111_poolB_beta" = Freq,
                  "p111_poolB_umi" = totalUMICount)

p111_betas_invitro %>% distinct(Beta_clonotype) %>% dplyr::count()
       n
1 288082

Convert to decimal frequency to %

p111_betas_pbmc <- p111_betas_pbmc %>%
  mutate(p111_pretreatment = p111_pretreatment * 100,
         p111_postvax = p111_postvax * 100,
         p111_prevax = p111_prevax * 100,
         p111_w40 = p111_w40 * 100)

p111_betas_invitro <- p111_betas_invitro %>%
  mutate(p111_DMSO_beta = p111_DMSO_beta * 100,
         p111_poolA_beta = p111_poolA_beta * 100,
         p111_poolB_beta = p111_poolB_beta * 100)

Join longitudinal and in vitro frequencies by beta chain

p111_betas <- full_join(p111_betas_pbmc, p111_betas_invitro)
Joining with `by = join_by(Beta_clonotype)`

Calculate fold changes across in vitro samples

p111_betas <- p111_betas %>%
  rowwise() %>%
  mutate(poolA_vs_dmso_fc = sum(c(p111_poolA_beta, pseudo_freq), na.rm = TRUE)/sum(c(p111_DMSO_beta, pseudo_freq), na.rm = TRUE),
         poolB_vs_dmso_fc = sum(c(p111_poolB_beta, pseudo_freq), na.rm = TRUE)/sum(c(p111_DMSO_beta, pseudo_freq), na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(log2fc_poolA_vs_dmso = log2(poolA_vs_dmso_fc),
         log2fc_poolB_vs_dmso = log2(poolB_vs_dmso_fc),
         max_log2fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ log2fc_poolA_vs_dmso,
                                             log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ log2fc_poolB_vs_dmso),
         max_fc_pool_vs_dmso = case_when(log2fc_poolA_vs_dmso > log2fc_poolB_vs_dmso ~ poolA_vs_dmso_fc,
                                         log2fc_poolB_vs_dmso >= log2fc_poolA_vs_dmso ~ poolB_vs_dmso_fc))

Save data

write.csv(p111_betas, "p111_betas_merged_Part1.csv", row.names = FALSE)

2.13 Get session info

sessionInfo()
R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Rocky Linux 8.10 (Green Obsidian)

Matrix products: default
BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so;  LAPACK version 3.9.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: America/New_York
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] lubridate_1.9.4 forcats_1.0.0   stringr_1.5.1   purrr_1.0.4    
 [5] readr_2.1.5     tidyr_1.3.1     tibble_3.2.1    ggplot2_3.5.1  
 [9] tidyverse_2.0.0 dplyr_1.1.4    

loaded via a namespace (and not attached):
 [1] gtable_0.3.6      jsonlite_1.8.9    compiler_4.3.2    tidyselect_1.2.1 
 [5] scales_1.3.0      fastmap_1.2.0     R6_2.6.1          generics_0.1.3   
 [9] knitr_1.49        htmlwidgets_1.6.4 munsell_0.5.1     pillar_1.10.1    
[13] tzdb_0.5.0        rlang_1.1.5       stringi_1.8.4     xfun_0.50        
[17] timechange_0.3.0  cli_3.6.3         withr_3.0.2       magrittr_2.0.3   
[21] digest_0.6.37     grid_4.3.2        rstudioapi_0.17.1 hms_1.1.3        
[25] lifecycle_1.0.4   vctrs_0.6.5       evaluate_1.0.1    glue_1.8.0       
[29] colorspace_2.1-1  rmarkdown_2.29    tools_4.3.2       pkgconfig_2.0.3  
[33] htmltools_0.5.8.1